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Abstract 



An algorithm to simulate the dynamics of a quantum state over a three-site lattice interacting with classical harmonic 
oscillators has been devised. 

The oscillators are linearly coupled to the quantum state and are acted upon by a fluctuation-dissipation process to 
take the equilibrium thermal environment into account, thus allowing to investigate how stochastic forces may affect the 
quantum dynamics. 

The implementation of the algorithm has been written in Ada95. 



lS The physical model 

Hosltein hamiltonian for a lattice connected to independent site oscillators can be distinguished into three terms 

h = h a + h b + h a b (1) 



o 
o 



corresponding to the primary system a, representing those states allowed to the quantum particle, the secondary system b, 
representing the vibrating lattice masses (m) and the interaction energy ab [?]. 

\£%n a three-site lattice, assuming the nearest neighbours approximation hold, the quantum states are described by the 
hamiltonian 

m : 

/i a = V(|l)<2| + |2>(l| + |2)<3| + |3>(2|) (2) 

wiiere \n) are the eigenstates of the isolated sites and V is the overlap integral between neighbour sites. 
■5^/ibrations are described by harmonic oscillators of equal mass m and frequency v, with respective position r n and 
cc^ugated variable p„, (where n indicates the location on the lattice): 

kb = i + m2cj2r i) + 2^ + + 2^ ( p z + m2uj2r D ( 3 ) 

InMraction terms are assumed to be linear with respect to r. 



h ab = e 1 |l)(l|+e 2 |2)(2|+e 3 |3)(3| (4) 



e^^= x r n , with x rea l positive, being the coupling energy at site n. Such physical model has been studued by Hennig [?] 

arfdiKenkre & Andersen [?], but without thermal environment. 

' 
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1.1 Lattice with free boundary conditions 

When loose ends boundary conditions apply to the primary system, we are lead to the following equation for the density 
matrix: 

Pii = iu (P2i - P12) 

P12 = -ILJ12P12 + iu (p 2 2 - P11 - P13) 

f>13 = -«Wi 3 pi3 + iuj (p 2 3 - P12) 

P21 = iui2P2i + iu a (pn + p 3 i - P22) 

f>22 = i^o (pi2 + P32 - P21 - P23) (5) 
P23 = -IU23P23 + i^o (Pl3 + P33 ~ P22) 

P3i = iuizpzi + iuj (p 2 i - P32) 

f>32 = i^23P32 + i^ (p22 - P31 - P33) 
p 33 = i(J (p 2 3 - P32} 

where lo — V/h is the frequency for the oscillation of the quasi-particle between two neighbour sites and u>kn = —w n k = 
(en - £fe) /fl = X (r n ~ r k ) /h. 
Using operators 



Ul 


= Pll - P22 


U2 


= P22 — P33 


Vi 


= i(Pl2-P2l) 


V2 


= i (P32 - P23} 


V3 


= i (P3i - P13) 


W\ 


= Pl2 + P21 


W2 


= P23 + P32 


W 3 


= Pl3 + P31 



(6) 



a new set of equations, that will be used to perform numerical simulations, can be obtained for the primary system. 

As a fluctuation-dissipation process has been finally attached to each oscillator, the complete dynamical model can be 
written: 

iii = -oj (2vi + v 2 ) 

u 2 = uj (vi + 2v 2 ) 

Vl = LU 12 Wl + lu (2ui + w 3 ) 

i>2 = -0J23W2 - u (2u 2 - w 3 ) 

v 3 = -W13W3 - io Q (wi - w 2 ) 

Wl = -W12V1 + UJ D V 3 

W 2 = U>23V2 ~ ^o«3 

W 3 = U13V3 - L0 o (vi + v 2 ) (7) 

fi = Pi/m 

r 2 = p 2 /m 

h = P3/m 

p\ = -mu?r x - - (c + u 2 + 2m) - 71P1 + fi (t) 

p 2 = -raw r 2 - - (c + u 2 - «i) -72P2 + h (t) 

p 3 = -raw r 3 - - (c - u-i - 2u 2 ) - 73P3 + h (t) 

where 7„ are the damping coefficients and /„ (t) models gaussian noise, with (5-shaped time correlation, satisfying the following 
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relations: 

having defined thermal energy 

(/cb being the Boltzmann constant); hence 

It can be readily verified that 



(/»(*)) = (8) 
(/»(*)/«(*')> = 27%_ t0 (9) 



k B T (10) 



(f n (t)) 1/2 = W) 1/2 (11) 



C = pll + P22 + P33 (12) 

is a constant of motion. A further constant, that will be used to monitor numerical simulations, is given by 

K = - (u\ + u 2 2 + Ul u 2 ) + v\ + uf + v\ + w\ + w\ + w\ (13) 

The dynamical problem is described by a system of stochastic real valued, first order, differential equations (0) that must 
be numerically solved. 

Notice that, restricting the number of sites to n = 1, 2 we are left with the spin-boson equation with noise: 





= —2lo vi 


Vl 


= U)i2W\ + 


Wi 


= -W12U1 


h 


= pi/m 


f-2 


= P2/m 


Pi 


= —mid r\ 


h 


= —muj 2 r2 



(14) 

|(c + «a + 2tti) + /i(t) 
- [C + U 2 - Ml) + h (*) 

that has been studied by many researchers; in particular, the author's work has been inspired the papers[?, ?], referenced to 
in the bibliography. 

2 Numerical integration method 

Stochastic differential equations of first order in time for a vector variabile x(£) take the from, 

x = F + gf (15) 

where F = (Fi) 1 = (-Fi(x, t)) 1 is the vector of deterministic fields, f = = (fi(t)) 1 the gaussian stochastic force, 

satisfying the conditions: 

(/(«)> - 

(/(t)/(f)> = 5(t-t>) (16) 

while th parameters set g = (gi) 1 n = {gi (x, t)) ln represents the interaction of the system with the thrmal bath. 

Equations given in ( |l5| ) are entirely coupled, since, in general, every components Fi and gi are functions of the whole 
variable x, 

The formal solution of (|l5|), assuming the implicit time dependence F = F(x(i)) e g = g(x(i)), is given by [?] 

Xi {t) - Xi(0) = f dt'Fi + f dt'gifi 1 = 1, ... n (17) 

Jo Jo 



3 



Given a time interval [0, h] that can be considered infinitely small for Fi e gi functions, their value at the arbitrary instant 
t' G [0, h], can be approximated by the first k terms of a Taylor expansion, centered at t = 0: 

Fi(i) = Ff + 5 K Fi (18) 
9i (t) = g° + 5*g t (19) 

The index k in ( |l^ , |l9|) indicates the truncation order of the series: 

5Fi = J2 l F *h ■ 5x i (0 + ^ E • fa i (O^* (0 + | E ■ ^ (*')<^ + ■ ■ ■ (20) 

% = E H • MO + 5 E Hfc • <M*')<^(*0 + ^ E • + • ■ ■ (2i) 

j " jk ' j'fc 

in which Fi and <?i, time derivatives evaluated at t = 0, have been represented by the following notation: 

Truncation at ft = the general solution obtained after substitution of ( |I^Jl9| ) into filq): 

Xi(h)- Xi (0)= [ dt'(F° + 5 K F l )+ [ dt' {g° i +5 K g i )f i (24) 
Jo Jo 



yields 



Xi (h) - Xi(0) = F°h + g° I dt'fi (25) 



where Ft = [Fi] = F° and gi — [tfe] — g° have been defined. 
Since the gaussian integral 



h 



Zi(h) = / dt'f(t') (26) 



is of h 1 ! 2 order, the second term in (25) is the lowest order approximation of the trajectory (|24|): 

fe| 1/2) = g?Z x (27) 

having defined 8xi = Xi (h) — Xj (0). 

Substituting ( p7| ) into ( pO| , ^T| ) , hence into ( |i~7| ) , and retaining only contribution of the order h in (|24|) , we obtain the first 
order correction 



5x 



(i) _ 



Jo j 

/•;'/' • E s*W I dt'z x {t')fi{t') 

f°H^[4. Sj % 2 (28) 



2 

3 
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The corresponding displacement on the trajectory is thus given by 

Sx.i — <5a4 2 + 8x\ 



(29) 



Again substituting the first order correction ( |2g| ) into the series ( 20 2lj) , hence into ( |l7| ) , then retaining only terms of the 
order ft, 3 / 2 in dgl), we get: 



5x 



(3/2) 



]T IF^ ■ 9°jZxh + [ dt> £ M 3 ■ SxVfitf) 

J2 m 3 ■ g^Zxh + £ dt' J2 ■ \F°h + [34 • 9lZl{t')h{t') 

[F^ ■ 9°jZxh + ]T [gilj ■ \F°h + \Y, [<?4 • 9% t dtZ 2 (t')Mt') 
j j V k J 

]T [ttlj • g]Z!h + ]T [<?4 • ( F°h + ± Y1 [ff4- • 9° k 'Z! ) 
j j \ ' k ) 



(30) 



In the simplified case in which g does not depend on x, the numerical algorithm to reach the order h 2 in approximating 
the exact solution of (FBI) can then be written 



t(h) - Xi (0) = g\Z x +F°h + J2 [Fi] s ■ 9]Z 2 + [F^ ■ [Fj] h 2 i = 1, . 



. . n 



or, in vectorial notation: 



x(/i) - x(0) = g°Zi + F°h + VFgZ 2 + -VF • Fh 2 



(31) 



(32) 



2.1 The algorithm for the three-site lattice 

Equations ([?]) belong to the class (|l5|); therefore the approximation method that has been recalled in the previous chapter 
can be implemented. 

To better expose our procedure, it is convenient switching to the vectorial notation; the system variable x, taking the 
form: 

Ul 
U-2 
i'l 
i' 2 
^3 
Wl 

w 2 

n 

r-2 

T3 
Pi 

P2 

_ P3 

evolves in time according to the equation 



x 



(33) 



x = F + G 



(34) 



5 



where F is the vector of the deterministic forces 



-cj (2vi + v 2 ) 
u> («i + 2v 2 ) 
uivzwi + uj (2ui + wa) 
-uj 23 w 2 - u (2u 2 - w 3 ) 

-WI3W3 - LU (lUl - W 2 ) 
—L0 12 Vl + LU a V 3 
L0 23 V 2 - ^o«3 

- L0 o (Vl + V 2 ) 

pi/m 
Pl/m 
Pz/m 



-mwVi - I (c + u 2 + 2ui) - 71P1 
-mui 2 r 2 — \ (c + u 2 



Ml) 



72P2 



3 

I (c - Ml - 2u 2 ) - 73P3 



while G contains the average values of the stochastic forces: 



G = 















(2 7 i0i) 

(2 7 20 2 ) 



1/2 
1/2 
,1/2 



(273C3 1 

The approximated solution up to the second order in time can be computed according to the formula (| 

x(/i) ~ x(0) + g • Zi + F • h + VF • g • Z 2 + ivF • F • h 2 
where Z\ and Z 2 are gaussian integrals and the force gradient is given by: 



VF = 









~2lo 


—UJ 
































uj 


2uj 


























2uo 














W12 





U! 




X2W1 











-2lu 














—^23 


u 





X2W 2 


-X3W2 




















—UJ 


U> 


-U13 


X1W3 





-X3W3 











-W12 
















X1V1 


-X2VI 

















UJ 23 
















-X1V2 


X3V2 













-CJ 


^13 











-X1V3 





X3V3 






































V m i 










































































-2xi/3 


-Xi/3 




















—m,\uj\ 








-71 


X2/3 


-X2/3 























— m 2 oj 2 








X3/3 


2X3/3 


























-m 3 o;| 
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The implementation of the algorithm according to the rule (|37|) is displayed in the last section . Ada95 language has been 
the language of choiche, since it is maintained that a rather easy and dependable parallel code, that could be required to 
perform long lasting simulations can be worked out in that environment, . 

Marsaglia & Tsang [?] fortran routine for random numbers generation has been used to compute Z\ and Z2 integrals. 

3 Results and Conclusions 

First run of the program allowed to plot variables : u\ = pn — P22 and u-i = P22 — P33, using the set of physical parameters: 

• inter-site interaction V = 0.09 

• free oscillators frequencies J 1 = 2a;10 -3 

• vibration-quantum state coupling \ — 4a;10~ 2 ; 

• 7 = 0.2 that imply lu 2 /2~/ 2 = 0.025 and x 2 /2mw 2 = 0.1 

Time integration step has been set to 10 -3 and the reliability of the computation has been monitored through the constant 
(]l3|), whose changes have kept below 10~ 7 . 




1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 




1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 

time 

(a) occupation differences ui = pn — P22 (lower) and 112 = P22 — P33 (upper) evolving in time 

Picture of the site occupation differences u\ = pn — P22 (lower) and U2 = P22 — P33 (upper) evolving in time (arbitrary 
units) shows the localized motion of the quantum state. 

The possibilty to model the dynamics of a three-site lattice coupled to stochastic oscillators allows one to study how 
a quantum state propagates in a real lattice, when vibrational motion can be described classically, then to investigate the 
duration of out of equilibrium states and the conditions under which energy transfer without dissipation may occur. 

4 Code listing 

Here included are the main tokens of code from the iterator procedure, in which the numerical integration at time step dh is 
performed: 
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procedure iterator (x : in out coordinates ; xh : out coordinates ; time_step : in double 

outer_cycle : in integer ; -- total number of trajectory -- points to print out 

v , ami , am2 , am3 , garni , gam2 , gam3 , oml , om2 , om3 : in double; chil , chi2 , chi3 , tl , t2 , t3 , elb , e2b , e3b : in dc 

output_file : in string ) 

is 



Fortran random numbers generation routine rnor is imported here: 

function rnor(dum : in Fortran. Fortran_Integer) return Fortran. double_precision; 
pragma Import (fortran, rnor, "rnor_") ; 



Definition of numerical parameters: 



dh : double := time_step ; -- time integratin step 

dh2 : double := time_step * time_step; -- squared time integration step 
dh25 := dh2 * 0.5; 



Definition of physical parameters: 

v2 : double := v * v ; -- 
akO := cl*cl + c2*c2 + c3*c3; 
qoml := oml*oml; 
qom2 := om2*om2; 
qom3 : = om3*om3 ; 
qgaml := gaml * garni; 
qgam2 := gam2 * gam2; 
qgam3 := gam3 * gam3; 

norm := (1.0/3.0)*( x(l)*x(l) + x(2)*x(2) + x(l)*x(2) ); 
norm := norm + x(3)*x(3) + x(4)*x(4) + x(5)*x(5); 
norm := norm + x(6)*x(6) + x(7)*x(7) + x(8)*x(8); 
norm := (3.0/4.0) * norm; 

Here the iteration starts: 

for i in 1 . . outer_cycle loop 
for j in 1 . . inner_cycle loop 



Definition of some useful expressionss: 



time 
wll 

¥12 

¥21 

¥22 

¥31 

¥32 

zll 

zl2 

z21 

z22 

z31 

z32 

oml 2 

om23 

oml 3 

qoml2 

qoml3 

qom23 

akey21 

akey23 

akey31 

ajayl 

ajay2 

ajay3 

omx3 : 

omxl : ; 



= time + time_step; 
double (rnor (dum)) ; 
double (rnor (dum)) ; 
double (rnor (dum)) ; 
double (rnor (dum)) ; 
double (rnor (dum)) ; 
double (rnor (dum)) ; 
all*¥ll ; 



a21*zll 
all*¥21 
a21*z21 
all*¥31 
a21*z31 
■■ del2 
de23 
■■ del3 
= oml 2 
= oml 3 
= om23 



+ a22*¥l2; 



+ a22*¥22; 



a22*¥32; 

- chi2*x(10) 

- chi3*x(ll) 

- chi3*x(ll) 

* oml2; 

* oml3; 

* om23; 

= - chi2*x(13)/am2 
= - chi2*x(13)/am2 
= - chi3*x(14)/am3 

= chil*(ak0+x(l)/2. 

= chi2*(ak0+x(l)/2. 

= chi3*(ak0-x(l)-x( 

■ om23 + oml3; 

■ oml2 + oml3; 



+ chil*x(9) ; 
+ chi2*x(10); 
+ chil*x(9) ; 



+ chil*x(12)/aml 
+ chi3*x(14)/am3 
+ chil*x(12)/aml 

0+x(2))/3.0; 

0-x(2)/2.0)/3.0; 

2)/2.0)/3.0; 
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Variables incrementation: 
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xh(l) := x(l) + 2.0*v*( 2.0*x(7) + x(6) ) * dh + 

v2*(x(2)-2.0*x(l)+3.0*x(5)) * dh2 + 
v*(2.0*om23*x(4) - oml2*x(3)) * dh2; 



xh(2) := x(2) - 2.0*v*( x(7) + 2.0*x(6) ) * dh + 

v2*(x(l)-2.0*x(2)-3.0*x(5)) * dh2 + 
v*(2.0*oml2*x(3) - om23*x(4)) * dh2; 



xh(3) := x(3) + ( oml2*x(6) + v*x(8) ) * dh - 

qoml2*x(3) * dh25 + v2* (x(4) -x(3) ) * dh25 + 
v*(oml2*x(2) + omxl*x(5)) * dh25 + 
akey21*x(6) * dh25 ; 



xh(4) := x(4) - ( om23*x(7) + v*x(8) ) * dh - 

qom23*x(4) * dh25 + v2* (x(3) -x(4) ) * dh25 + 
v*(om23*x(l) - omx3*x(5)) * dh25 + 
akey23*x(7) * dh25 ; 



xh(5) := x(5) - ( oml3*x(8) + v*x(6) + v*x(7) ) * dh - 

qoml3*x(5) * dh25 + v2*(x(l)-x(2) - 2.0*x(5)) * dh25 + 
v*(omxl*x(3) - omx3*x(4)) * dh25 - 
akey31*x(8) * dh25 ; 



xh(6) := x(6) - ( oml2*x(3) - v*x(2) - v*x(5) ) * dh - 

qoml2*x(6) * dh25 - v2* (5 . 0*x(6)+3 . 0*x(7) ) * dh25 - 
v*omxl*x(8) * dh25 - 
akey21*x(3) * dh25 ; 



xh(7) := x(7) + ( om23*x(4) - v*x(l) + v*x(5) ) * dh - 

qom23*x(7) * dh25 - v2* (3 . 0*x(6)+5 . 0*x(7) ) * dh25 - 
v*omx3*x(8) * dh25 - 
akey23*x(4) * dh25 ; 



xh(8) := x(8) + ( oml3*x(5) - v*x(3) + v*x(4) ) * dh - 

qoml3*x(8) * dh25 - 2.0*v2*x(8) * dh25 - 
v*(omxl*x(6) + omx3*x(7)) * dh25 + 
akey31*x(5) * dh25 ; 



xh(9) := x(9) + (x(12)/aml) * dh - 
qoml*x(9) * dh25 - 
dh25 * (ajayl + gaml*x(12) )/aml + 
(Dl/aml) * zl2 ; 



xh(10) := x(10) + (x(13)/am2) * dh - 

qom2*x(10) * dh25 - 

dh25 * (ajay2 + gam2*x(13))/am2 + 

(D2/am2) * z22 ; 



xh(ll) := x(ll) + (x(14)/am3) * dh - 

qom3*x(ll) * dh25 - 

dh25 * (ajay3 + gam3*x(14))/am3 + 

(D3/am3) * z32 ; 



xh(12) := x(12) + Dl*zll - 

( aml*qoml*x(9) + ajayl + gaml*x(12) ) * dh + 
qoml*(aml*gaml*x(9) - x(12)) * dh25 + 
(chil*v*x(6) + gaml*ajayl + qgaml*x(12)) * dh25 - 
gaml*Dl*zl2; 



xh(13) := x(13) + D2*z21 - 

( am2*qom2*x(10) + ajay2 + gam2*x(13) )*dh + 
qom2*(am2*gam2*x(10) - x(13)) * dh25 + 

(-chi2*v*(x(7)+x(6)) + gam2*ajay2 + qgam2*x(13)) * dh25 - 
gam2*D2*z22; 



xh(14) := x(14) + D3*z31 - 

( am3*qom3*x(ll) + ajay3 + gam3*x(14) )*dh + 
qom3*(am3*gam3*x(ll) - x(14)) * dh25 + 
(chi3*v*x(7) + gam3*ajay3 + qgam3*x.C14) ) * dh25 - 
gam3*D3*z32; ^ 



for ii in 1 . . conf ig_space loop 

x(ii) := xh(ii) ; 
end loop; 



end loop; 



end loop; 
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